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I. INTRODUCTION 



The invisible axion is one of the best motivated candidates for cosmic dark matter. 
The axion is the pseudo-Nambu-Goldstone boson resulting from the spontaneous break- 
ing of a U(l) global symmetry known as the Peccei-Quinn, or PQ, symmetry. The PQ 
symmetry is introduced to explain the apparent smallness of strong CP-violation in QCD 
l]J. Although there are other possible solutions to the strong CP problem 0, and the 
origin of the axion in the breaking of a global symmetry has been criticized the axion 
remains the best known cure for the disease of strong-CP violation. 

There are stringent astrophysical and cosmological constraints on the prop- 
erties of the axion. In particular, the combination of cosmological and astrophysical 
considerations restrict the axion decay constant f a and the axion mass m a to be in the 
narrow windows 10 10 GeV < f a < 10 12 GeV, and 10" 5 eV < m a < KT 3 eV §. The 
contribution to the mean density of the Universe from axions with mass in this win- 
dow is guaranteed to be cosmologically significant. Thus, if axions exist, they will be 
dynamically important in the present evolution of the Universe. 

In addition to the usual role in the evolution of primordial density fluctuations and 
the formation of large-scale structure common to all cold dark matter candidates, ax- 
ions have unique features as dark matter. The energy density in axions corresponds 
to coherent scalar field oscillations, driven by a displacement of the initial value of the 
field (the "misalignment" angle) away from the eventual minimum of the temperature- 
dependent potential. During the QCD epoch fluctuations in the misalignment angle on 
scales comparable to the Hubble radius at that time || are transformed into large ampli- 
tude density fluctuations, which later lead to tiny gravitationally bound "miniclusters" 
||. It was found that the density in miniclusters exceeds by ten orders of magnitude the 
local dark matter density in the Solar neighborhood This might have a number of 
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astrophysical consequences, as well as implications for laboratory axion searches JO 

In previous studies of the evolution of the axion field around the QCD epoch, the 
effect of spatial gradients of the axion field were either neglected, or were included in a 
limit where the non-linear potential was approximated by a linear harmonic potential. 
Both approximations are adequate for temperatures well below the QCD scale where 
the coherent axion oscillations can be treated as pressureless, cold dust. However, in a 



previous paper |LI[ we found that just at the crucial time when the inverse mass of the 
axion is approximately the size of the Hubble radius and fluctuations of misalignment 
angle are still of order tt, both the non-linear interaction and the gradient terms are 
important, and a full field-theoretical calculation is needed. Here, we present the results 
of a 3-dimensional numerical study of the evolution of the inhomogeneous axion field 
around the QCD epoch. We find that the resulting axion clumps are much denser than 
previously thought, even reaching the critical conditions for Bose star formation jEJ. 

In Sec. II we review the basic scenario for the evolution of the axion field around 
the QCD epoch. In Sec. Ill, after deriving the equations of motions in a suitable form, 
we present the results of (3 + l)-dimensional numerical calculations of an initial white- 
noise axion distribution. We find that the non-linear potential leads to the formation 
of dense, roughly spherical, soliton-like axion configurations we call axitons. We then 
follow the subsequent evolution of these spherically symmetric configurations in a (1 + 1)- 
dimensional calculation. Sec. IV is devoted to the consideration of an initial axion field 
that results in a network of topological defects. We discuss how the usual picture of 
axion strings slicing up axion domain walls is modified by the inclusion of the non- 
linearities of the true axion potential. We find that rather than the simple picture 
of axion strings destroying walls by punching holes in them, unstable pseudo-breather 
solitons are formed which decay to axitons. In the final section we discuss some possible 
physical consequences of very dense axion clumps. 



II. COSMOLOGICAL EVOLUTION OF THE AXION FIELD 



The axion story begins with PQ symmetry breaking. This symmetry breaking occurs 
when a complex scalar field with non-zero PQ charge develops a vacuum expectation 
value. This PQ symmetry breaking can be modeled by considering a potential of the 
standard form V($) = A (\(f\ 2 — f 2 /2j . The axion is the resulting Nambu-Goldstone 
degree of freedom from spontaneous breaking of the global symmetry. After PQ sym- 
metry breaking at T ~ f a , but before QCD effects are important, the axion is massless. 
However since the PQ symmetry is anomalous, it is broken explicitly by QCD instanton 
effects, leading to a mass for the axion. In general the instanton effects respect a residual 
Zn symmetry, and the axion develops a potential due to instanton effects of the form 
V(a) = m 2 (f a /N) 2 [l — cos(Na/ f a )}. The axion field is often represented in terms of an 
angular variable 9 = Na/f a , and if 9 is taken as the dynamical variable, its potential for 
AT = 1 is 

V{9) = m 2 a (T)f 2 (l - cos9) = A*(T)(1 - cos£). (2.1) 

Because QCD instantons are large, with a size set by Aqq D , their effects are strongly 
suppressed at high temperatures. So for T ^> Aqcd, the axions are effectively massless. 
For T ^> Aq CD , the temperature dependence of the axion mass scales as 

m 2 a (T) = m 2 a (T*)(T/T*)- n ; n = 7.4 ± 0.2. (2.2) 

When the field 9{x) is created during the Peccei-Quinn symmetry breaking phase 
transition at T ~ f a , it should be uncorrelated on scales larger than the Hubble radius 
at that time [0. As the temperature decreases and the Hubble radius grows (in a 
radiation-dominated Universe the Hubble radius grows as Rh{T) = H^iT) oc T~ 2 ), the 
field becomes smooth on scales up to the Hubble radius. This continues until T = T± ~ 1 
GeV when the axion mass "switches on," i.e., when m a (Ti) w 3H(Ti), and the axion 
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mass begins to become important in the equations of motion. Coherent axion oscillations 
then transform fluctuations in the initial amplitude of the axion field into fluctuations in 
the axion density. 

Since the initial amplitude of the coherent axion oscillations on the scale of the Hub- 
ble radius is uncorrelated, one expects that typical positive density fluctuations on this 
scale will satisfy p a ~ 2p a , where p a is mean cosmological density of axions @. At 



the temperature of equal matter and radiation energy density, Teq = 5.5Q a h eV [|14 
non-linear fluctuations will separate out as miniclusters with pmc ~ lCr 14 gcm -3 ||. 
The minicluster mass will be of the order of the axion mass within the Hubble radius 
at temperature 7\, Mmc ~ 1O~ 9 M . The radius of the cluster is -Rmc ~ 10 13 cm, and 
the gravitational binding energy will result in an escape velocity of v e /c ~ 10~ 8 . Note 
that the mean phase-space density of axions in such a gravitational well is enormous: 
n ~ Pa m- 4 v; 3 ~ 10 48 /i 4 2 , where / 12 = f a /10 12 GeV. 

We will show below that due to non-linear effects, a substantial number of regions 
at Aqcd > T > Teq can have an axion density orders of magnitude larger than 2 p a . 
These form because the non-linear effects in the axion potential lead to the formation of 
pseudo-soliton objects we call axitons. 

The axitons are not true solitons because the field coherently oscillates inside the 
axiton. The oscillations of the field lead to a red-shift of the energy density of the field 
in an expanding Universe. Quantitatively, axitons resemble breathers of the (1 + 1)- 
dimensional sine-Gordon model. 

Eventually the energy density of the axiton is red-shifted to sufficiently small values 
of the axion field so that non- liner ait ies can be neglected, and the axiton configuration 
is frozen in the cosmological expansion as is any linear fluctuation. However the energy 
contrast relative to the homogeneous background will be large. 

Once an axiton forms, its energy density scales as T 3 for T > Teq, so we can write 



Paxiton = 3 (1 + $)Teqs/4, where $ depends upon the initial conditions of the axion 
field, i.e., the misalignment angle and its gradients at T\. Here, s is the entropy density, 
and $ = corresponds to the mean axion density. The energy density inside a given 
fluctuation is equal to the radiation energy density at T = (1 + $)T E q. At that time the 
self gravity of the fluctuation comes to dominate, and if $ > 1 it separates out from the 
cosmological expansion, collapses, and forms a minicluster with density^] 



Pmc - 140$ 3 (1 + $)p a (T EQ ) « 3 x l(T 14 $ 3 (i + $) (tt a h 2 ) g cm" 3 , (2.3) 
Even a relatively small increase in $ is important because the final density depends upon 

Ours is not the first proposal that non-linear effects can lead to large values of $. 
One mechanism whereby non-linear effects can lead to amplification of the axion density 
was recognized in Ref . || . In that analysis it was proposed that some correlation regions 
can have values of <3> larger than one because the closer the initial value of 9 is to the top 
of the axion potential, the later axion oscillations commence. However, this effect alone 
is not very significant. If the closeness of the initial angle to the top of the potential 
is parameterized by £ = (tt — O^/n, then for £ in the range 0.1 < £ < 1CT 3 , $ — 1 ~ 
1.5(#j/7r) 2 £ -a35 , and $ is significantly larger than 1 only for initial values very finely 
tuned to the top of the potential |Tl|. Moreover, the axion field is not exactly coherent 



on scales of the Hubble radius, and even small fluctuations will spoil this simple picture. 

Our scenario for the generation of axion miniclusters mainly depends upon the in- 
terplay of the non-linear effects in the potential and gradients in the axion field. The 
interplay of these two effects will lead to the formation of axiton configurations in the 
axion field. At temperatures T 3> T\, the potential is negligible in the equations of mo- 
tion compared to the gradient terms which force the field to be homogeneous on scales 



The factor of 140 results from a detailed calculation. 
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less than the Hubble radius. At T <C Aq C d, gradients can be neglected and one can 
treat the evolution of fluctuations as that of a pressureless gas. Clearly, around the QCD 
epoch when the potential just starts to become important in the equations of motion 
the gradient terms are still important, and since the initial amplitude can be close to 
7r, the non-linear nature of the potential is also important. In order to find the energy 
density profile at freeze out one has to trace the non-linear inhomogeneous field evolution 
through the epoch T ~ 7\. 



III. INHOMOGENEOUS AXION FIELD EVOLUTION 

A. Equations of Motion 

We start with deriving the equations of motion for the axion field in a form suitable 
for numerical calculations. In an expanding, spatially flat Universe with scale factor R(t), 
the equation of motion for the axion field takes the familiar form 

6 + 3 f ^ ~ 1W) A26 + m ' (t) sin ^ = °' (3,1) 

where dot denotes time derivative and A is the Laplacian with respect to comoving 
coordinates x. 

Rather than cosmological time, it is convenient to work with a conformal-time co- 
ordinate. In a radiation-dominated Universe the conformal time is proportional to the 
scale factor R. Using R as the independent variable, the equation of motion is 
d 2 9 2 d6 1 1 , 2 1 2 



AO + ^m 2 a (R) sin 9 = 0. (3.2) 



dR 2 RdR R?R?{t) R* 
Using the Friedmann equation, along with the dependence of the expansion rate upon 
R in a radiation-dominated Universe, we can express R 2 in terms of the Hubble radius 
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at some arbitrary epoch (denoted by subscript 1): R 2 = H 2 R 2 = H 2 (Ri)Rf/R 2 . Now 
defining conformal time 77 as r\ = R/Ri, the equation of motion is 
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where prime denotes d/dr). We use Eq. ( |2.2| ) to find that in conformal time the mass 
evolves as m 2 (R) = m 2 (Ri)i] n . We can use the remaining freedom in the choice of R± 
to simplify the equation of motion by making the choice m 2 (Ri) = H 2 (Ri), i.e., rj = 1 
corresponds to the epoch when the inverse of the axion mass is equal to the Hubble 
radius. The equation of motion then takes the form 

6" + -9' - A 2 9 + r] n+2 sin 9 = 0, (3.4) 
V 

where A is now the Laplacian with respect to comoving coordinates, x = H(Ri)Rix. 
In other words, x = 1 corresponds to the Hubble radius at the epoch when the Hubble 
radius is the inverse of the axion mass. 

The equation of motion can be written as a wave equation by the introduction of the 
field ip =. rj9: 

•ij) - A 2 tp + r] n+3 sm(i/;/r]) = . (3.5) 

The equation of motion is finally in a form convenient for the study of the evolution 
of the axion field during the epoch when the mass switches on. In Table I we give the 
scaling with rj of several important physical length and mass scales. We next turn to the 
specification of the initial conditions. 

B. Initial Conditions 

At rj <^ 1 the potential term in Eq. (|3.5|) can be neglected, and the solution of 
the wave equation can be expressed simply as a sum of Fourier harmonics. As usual, 
there will be two sums over frequency u: one sum proportional to sin(u;?7) and one sum 



Table I: The scaling of physical quantities with conformal time r\. To find the scaling in coordinate 
distance, a length must be divided by 77, and a mass multiplied by rj. 



TIME 


t(v) = 


t(v = 


l)rf 




TEMPERATURE 


T( V ) = 








SCALE FACTOR 


R(v) = 


--R(v 


= l)r? 




AXION MASS 


m a {rj) 


= m a 


(77 = l)?7 n / 2 n 


= 7.4 ±0.2 


HUBBLE RADIUS 


Rh(v) 


= H 


- 1 ( V )=R h (t 1 = 





proportional to cos(cj?]). In the decomposition of the 6 field, terms like A{uj) sm.{urj) / (cur}) 
and B(uj) cos(curj)/(cur]) will appear. Assuming a finite amplitude for fluctuations of 
6 (of order of several it) on scales larger than the Hubble radius at the epoch of the 
Peccei-Quinn phase transition, we see that the coefficients B(cu) must be proportional 
to Aqcd//<i, while the coefficients A(cu) are of order unity. In other words the terms 
proportional to cos(u;y7) correspond to decaying modes on scales larger than the Hubble 
radius and can be neglected. Finally, assuming that on large scales the distribution for 
6 is white noise, we obtain 

9 = An sin(piX + ip lijk ) smfay + (p 2ij k) sin(p k z + (p 3ijk ) , (3.6) 

ijk 

where </?'s are random phases and ou 2 = pf + + p\. On scales larger than the Hubble 
radius the field distribution is frozen, while modes smaller than the Hubble radius are 
redshifted away. 

We numerically evolved this distribution starting from initial time 77 = 0.4 in a box 
of size L = 4 with periodic boundary conditions. There were 100 3 grid points in the box. 
Each of the momenta in the field decomposition took six discrete values, p n = 2im/L, 
with n — 1, ...,6. So, in total there were 3 x 6 3 random phases, each with values in the 
interval < (p < 2ir. 

The final parameter to be chosen is the magnitude of A. Recall that for N = 1, the 
axion potential is periodic with period 2n. We will consider two possibilities: A = 1 and 
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A = 2. For the case A — 1 it is unlikely that domain walls will form in a box of the 
size we study. However for A = 2 (the more physical choice) domain walls are produced 
at about 1 per horizon. We will present some results where A = 2, but for the most 
part we will consider in detail calculations with the A = 1 initial condition, since we are 
interested in the structure of density enhancements that are not associated with axion 
domain walls. So unless otherwise specified, our results will be for A = 1. 

The initial conditions are illustrated in Fig. 1 by a 2-dimensional slice through the 
3-dimensional box. The height above the plane is proportional to the axion energy 
density. Since at this epoch the axions are relativistic, it is convenient to scale their 
energy density by rf. The energy density shown in Fig. 1 is scaled by rf / p a (rj = 3) where 
p a {v) represents mean axion energy density at a given r\. Note that the Hubble radius at 
this epoch [r\ = 0.4) is 0.4 in the units of the figure, and the inverse of the axion mass is 
75 units. 

C. Results of Numerical Calculations 

1. (3 + 1)- dimensional evolution 
We first present the results of numerical calculations with A = 1, where density peaks 
arising from collapsing domain walls are filtered out so as to isolate the effects due to 
axitons. In order to present the results of the calculations we will take a two-dimensional 
slice through the three-dimensional box, and plot the energy density as the height above 
the plane. We have analyzed the time evolution of the energy density in several different 
slices. All of the slices generally look alike. The most important (and generic) feature is 
the development of large-amplitude peaks. As the system evolves in time, the peaks in 
the energy density, the axitons, increase in magnitude and become more compact. We 
present the results in the z = const plane, which intersects the point with the maximum 
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energy density. We emphasize that all slices through the box are quantitatively similar. 
We normalize the energy-density by comparing it to the energy density of a homogeneous 
axion field, p a {v)i with initial amplitude equal to the rms value of the misalignment angle, 

Orms =71"/ \/3. 

In order to isolate the effect of the non-linearities in the axion potential, we also evolve 
the same initial conditions with a harmonic axion potential, V(6) = m 2 (T)f 2 6 2 /2, and 
compare the evolution of the harmonic potential model to the axion model. 

The distribution of the axion energy density in the reference plane at time corre- 
sponding to 77 = 2 is shown in Fig. 2a for the harmonic potential, and in Fig. 2b for 
the axion potential. The maximum energy density peak that picks the reference plane is 
clearly seen in Fig. 2b, its top portion is choped off to fit overall the scale of the figure. 

The distribution of the axion energy density in the reference plane at time corre- 
sponding to i] = 3 is shown in Fig. 3a for the harmonic potential, and in Fig. 3b for 
the axion potential. Again, the tops of the four peaks in Fig. 3b are chopped off; their 
heights are in excess of 100! Of course the peaks are only evident for the axion potential 
model. 

Comparing Fig. 3b to Fig. 2b, we see that for the axion potential most of the high 
magnitude peaks grow considerably in height and became thinner, while most of the 
low amplitude peaks remain almost unchanged, i.e., they are in the linear regime and 
consequently are frozen by the cosmological expansion. There are some peaks (some even 
relatively high at 77 = 2) which decreased in amplitude. Those peaks represent the tales 
of the density clumps which reach their maximum at some other value of z. All high 
density peaks contract in the coordinate volume, those which decreased in height simply 
moved out of our reference plane. High density peaks do not develop in the evolution of 
the harmonic potential, and the evolution proceeds as was assumed in the linear analysis 
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There is insufficient resolution on this grid to proceed further in time with the axion 
potential, but the harmonic potential can be evolved further. In Fig. 4 we present the 
result of the distribution of the axion energy density in the reference plane at time 
corresponding to rj = 4 for the harmonic potential to demonstrate that as expected the 
evolution of the field in the linear regime is frozen by 77 = 3. Note that the typical 
magnitude of the peaks is about 2 for the harmonic potential. 

There is a simple, heuristic explanation for the fact that non-linear effects lead to 
the formation of high density peaks. The average pressure over a period of homogeneous 
axion oscillations in the axion potential potential is negative,^ and is equal to (P) ~ 
— A^(T)6'q/64, where 80 is the amplitude of the oscillations |16| (this formula is valid for 
9 <C n ; as #0 — > 7T, the field spends more and more time near the top of the potential, 
and (P) — > — 2A^). In other words, the axion self- interaction is attractive. The larger 
the amplitude of oscillations inside the fluctuation, the more negative the pressure inside, 
and consequently, fluctuations with excess axions will contract in the comoving volume. 
In addition, matter with a smaller pressure suffers less redshift in cosmological expansion. 

Before continuing our exploration of the evolution of the peaks by means of a 1- 
dimensional calculation, we present some results of calculations with A = 2, where do- 
main walls are much more likely to form than the above calculation with A — 1. The 
best way to illustrate the presence of domain walls is by a contour graph, where the 
shading represents the amplitude of the axion energy density. We show a graph of the 
energy density distributions for the axion potential at time rj = 2 with A = 2 in Fig. 5a 
and compare it to a similar contour graph for A = 1 at rj = 3 in Fig. 5b. In Fig. 5a two 
shells of collapsing domain walls are clearly visible in the lower left hand corner and in 

the upper right hand corner. Such configurations do not appear in Fig. 5b. The density 

2 Of course the average pressure is dominated by relativistic species at this time. It is the pressure 
contributed by the axions that is negative. 
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peaks in Fig. 3b, the axitons, are not related to axion domain walls. 

In order to learn the fate of the high density peaks, we have choosen one of them 
in Fig. 3 and generated the corresponding spherically symmetric initial conditions at 
i] = OA and evolved it in time. We now describe the result of this calculation. 

2. (1 + 1)- dimensional evolution 

The axiton we choose to examine is the one near the center of the grid of Fig. 3, with 
grid coordinates {2.24, 1.92} (the grid coordinate of the plane of Fig. 3 in the z-direction 
is 1 . 76 — also near the center of the 3-dimensional box — and the axiton we chose is almost 
at its maximum in this plane, having an absolute maximum at z = 1.80). 

The dependence of the axion field upon time at the reference point at the center of 
this peak in our (3 + l)-dimensional numerical calculation is shown in Fig. 6 by the solid 
curve. We can compare this evolution to the evolution of a massless axion field with the 
same initial conditions since we are able to calculate its evolution anaytically from the 
massless wave equation with initial conditions given by Eq. (|3.6|). The evolution of a 
massless field at the reference point is shown in the Fig. 6 by the dashed line. 

It is then straightforward to construct a spherically symmetric solution to the mass- 
less counterpart of Eq. (|3.5|) which has exactly the same time dependence as Fig. 6 in the 
center of the configuration. We start with the field decomposition, Eq. (|3.6| ), substitute 
the values of the coordinates of the given spatial point, and multiply each time harmonic 
by sin(u;r)/W. So the dashed line in Fig. 6 also represents the time dependence for a 
massless field at the reference point of the (3 + l)-dimensional calculation and also at 
the center of a (1 + l)-dimensional spherically symmetric configuration. (Away from the 
reference point even the massless field will evolve differently in the (l + l)-dimensional cal- 
culation and the (3 + l)-dimensional calculation.) We can use the resulting configuration 
9 = 9(j], r) for generating spherically symmetric initial conditions which will approximate 
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the peak of our choice for the runs with the actual axion potential. The result of such a 
calculation with A = 1 is shown in Fig. 6 by the dotted line. 

This plot is instructive in the evaluation of the accuracy of the numerical code. At 
i] < 1 the axion field is approximately massless, and the dotted curve is indistingushable 
from the dashed curve (there were 10 4 grid points in r-direction in the case of spherical 
symmetry) and the solid curve deviates only very slightly (near the extrema) from the 
dashed curve. This suggests that the use of 100 grid points in each spatial direction in 
the (3 + l)-dimensional calculations is adequate. 

From Fig. 6, we see that the axion mass effectively switches on at a time r] ~ 1.3. By 
this time the amplitude of the massless field is greater than unity.0 This means that for 
the evolution of the field using the axion potential the oscillations start in the non-linear 
regime 9 > 1 in the region that will develop into a axiton. We see also that the non- 
linearity is strong enough to force the density peak to collapse not only in coordinate 
space, but also in physical space as well, since by the time r] ~ 3 the amplitude of axion 
field oscillations are in the non-linear regime and growing (this is somewhat difficult to 
see in the figure). The rate of growth in the non-spherical case is much slower (compare 
the solid and dotted lines). This makes sense because we expect a spherical collapse to 
lead to a denser central region. 

So in general, there are two competing effects in the evolution of the axion field. 
1) A contraction of the axiton due to the pressure difference, leading to an increase in 
amplitude in the center, and 2) a decrease in amplitude of the oscillations due to the 
expansion of the Universe. We found that with a sufficiently large initial amplitude at 
the start of oscillations, the first process wins, and the amplitude in the center of the 

axiton increases to values 9 > n. In the opposite case, e.g., if the initial amplitude when 

3 Note that around rj — 2 the amplitude for the massless calculation is slightly larger than the ampli- 
tude of the calculations using the true axion potential. This is because as the axion mass switches on 
the amplitude of the axion field decreases. 
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oscillations commence is in the linear regime, the amplitude monotonically decreases in 
time. 

In realistic axion models the axion mass does not continue to grow with rj, but satu- 
rates at its zero-temperature value around T ~ Aqcd- For axion masses in the allowed 
window this corresponds to rj > 6. Due to the steep power-law dependence of the func- 
tion 777.(77) m Eq. ( |2.2|) , the period of the field oscillations becomes very small by 77 ~ 6, 
and direct numerical methods fail even in the case of spherical symmetry. In order to 
follow the evolution of the fluctuation up to freeze out we must assume that the mass 
saturates at a smaller value of 77. This would correspond to a larger value of f a . We have 
approximated the procces of axion mass saturation by the the simple formula 



taking 7/ c = 3.5. This value of t/ c corresponds to too large value of f a ~ 4x 10 13 GeV, which 
would give Vt a h 2 in excess of one. However we expect that qualitatively the evolution of 
the axion field will have the same basic properties for larger values of 77 c (smaller values 



We can vary the initial overall amplitude A of our spherically symmetric configu- 
rations. This has the effect of spanning different initial conditions of a well defined 
one-parameter family of axitons. Moreover, varing A is easier than choosing different 
peaks in Fig. 3. 

The time dependence of the field in the center of the fluctuation that will develop 
into an axiton for A = 0.73 is presented in Fig. 7, and for A = 0.77 in Fig. 8. In 
both cases the configuration collapses and the amplitude of 9 rapidly increases in the 
center, even exceeding the value of ir. This is followed by a period of several rebounds. 
An expanded view of the rebounds is shown in Fig. 9. During each rebound (eight in 
total in both cases) relativistic axions are emitted. We can see the signature of axion 




(3.7) 



Of fa). 
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emission by looking at the radial profile of an axiton. In Fig. 10 we show the profile of 
the axiton of Fig. 9 at three instants in time during one oscillation period. The emission 
of relativistic axions is seen in the outgoing waves of Fig. 10. The emission of relativistic 
axions reduces the energy of the central configuration below some critical value, at which 
point a pseudo-soliton, an axiton, is produced. 

The radial coordinate r is the spherical analogy to x. At the values of r\ in Fig. 10 
the axion mass has saturated to m a = 3.5 3 - 7 r] ~ 100?y. Therefore in the units of Fig. 10, 
the Compton wavelength of the axion is O.OI77 -1 , and the axiton is obviously much larger 
than m" 1 — it is indeed a soliton-like configuration. 

The axiton is a quasi-stable (on time scale m~ l ) solution of the field equations in an 
expanding Universe. Since there are no absolutely stable spherically symmetric breather- 
like solitons in flat space, in Minkowski space-time an axiton configuration will gradually 
decay anyway without the emission of axions present in the violent oscillations seen in 
Figs. 7 and 8. In an expanding Universe the situation is different. Once the axiton enters 
the linear regime it becomes frozen by the cosmological expansion, and behaves as a 
clump of coherent field oscillations (or ultra-cold axions). 

The final energy density profile of this configuration for the case A = 0.77 is shown 
in Fig. 11. At time 77 = 9 (dotted line), outgoing secondary waves are still seen in the 
tail of the configuration. By time 77 = 11 there is no evidence of outgoing radiation. The 
amplitude of the energy density at r = is 23.5 at 77 = 9 and 12.9 at rj = 11 (the energy 
density in this graph is not normalized to the homogeneous background). It is clear that 
the energy density in the center scales as ?7~ 3 (e.g., 23.5/12.9 = (11/9) 3 ), confirming that 
the linear regime has been reached, the fluctuation is frozen, and the number of axions 
per comoving volume is conserved. The energy density of a homogenoeus background at 
i] = 10 with i] c = 3.5 and initial amplitude equal to the rms value of 6 is 0.85 in the units 
of the figure. Thus, the fluctuation of Fig. 11 has an energy density contrast of 20. 
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Not all fluctuations that pass through the non-linear regime contract in physical 
space. For example, a sample spherical fluctuation for A = 0.70 does not collapse. 
The corresponding energy density profiles of this fluctuation at two moments of time 
are presented in Fig. 12 by the dashed lines. This should be compared to the solid 
lines, which are the energy density profiles for a fluctuation with A = 0.73 which does 
undergo collapse. We see that the slope of the energy density in the non-linear tail 
tends to a power law p oc r~ 3 prior to the collapse. This leads to an increase in field 
amplitude in the center, while due to the overall expansion of the Universe, the amplitude 
decreases. For A = 0.73, the first process wins for some period of time, see Fig. 7, while for 
A = 0.7 the general expansion dominates, and the amplitude of the oscillations decreases 
monotonically. However, the decrease in amplitude is much slower than it would be with 
the harmonic potential, and the final energy density contrast with rj c = 3.5 and A = 0.7 
is 45. 

For comparison we also present in Fig. 12 the energy density profile of the fluctuation 
with A = 0.77 at rj = 11. Remarkably, it has the same power-law slope, p oc r~ 3 , despite 
the fact that this profile represents a fluctuation that has undergone "violent oscillations" 
accompanied by axion emission (see Figs. 8, 10 and 11). 

Since the axion interaction is attractive, one can expect that bound states of axions 
can form. One example of such a bound state is the well known "breather" solution in 
the (1 + l)-dimensional sine-Gordon model. In (3 + 1) dimensions this solution possesses 
planar symmetry and turns out to be unstable with respect to fragmentation (we dicuss 
this further in the next section). If a spherically symmetric counterpart of the "breather" 
would exist in Minkowski space-time, it would behave in an expanding Universe just as 
the fluctuation shown in Fig. 11. Thus the axiton is related to a spherically symmetric 
breather. 

Suppose we can extrapolate these results to the range of realistic axion models, i.e., 

16 



to larger values of rj c corresponding to smaller values of f a . Then we must consider 
the possibility of producing enormous density contrasts. Indeed, both the increase in 
axion mass and the expansion of the Universe adiabatically decreases the amplitude of 
axion oscillations in the linear regime (or in the homogeneous state), so that at T < 
Aqcd the corresponding background energy density is about p a ~ TeqT 3 , where Teq ~ 
5.5Q a h 2 eV is the temperature of equal radiation and axion energy density In the case of 
a collapsing non-linear fluctuation, the final field configuration is the output of non-linear 
dynamics. Let ^ be the amplitude of field oscillations in the axiton at the time when 
it enters the linear regime at Tl = Tx/rji. Then the corresponding energy density in the 
fluctuation will be at this time about A 4 #|. The ratio of the axiton energy density to 
the homogeneous background axion energy density will be 

1 + $ « A A a el V l/T EQ Tl (3.8) 

Using the results from Figs. 7 and 8 (Ol ~ 0.1, and r\ L > 6),Q we obtain 1 + $ « 10 4 prior 
to gravitational decoupling of the fluctuations from the cosmological expansion. 

Although this possibility is exciting, a word of caution is necessary. Non-linear dy- 
namics is rather unpredictable, and one can not exclude the possibility that at r] c > 6 
all collapsing non-linear fluctuations somehow dissipate, leaving very small 6l- Note also 
that non-spherical configurations can evolve quite differently than the spherical configu- 
rations. 

The range of initial conditions which will lead to monotonic behavior of the amplitude 
in the non-spherical case is expected to be wider. Our point of view is that spectrum of 
energy density contrasts can span the entire range from order 1 up to of order 10 4 or even 
larger. However, at this time we have nothing to say in regard to the number density of 
peaks function of its amplitude. 



4 



In any case, rj^ will be larger than r/ ci the value of r\ where the axion mass saturates to its zero- 



temperature value [see Eq. (3.7)], and r/L > 6 seems a very conservative estimate. 
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So far we have neglected the presence of other non-linear structures which can be 
formed by the axion field during the QCD epoch, namely axion domain walls and walls 
bounded by strings. We now turn to the question of their fate and their contribution to 
the dark matter distribution. 



IV. AXION BREATHERS 

In general, there are four sources of cosmic axions. The first source is thermal axions 
18fl . The second source, related to the initial misalignment of the axion degree of freedom 



from its true minimum, was discussed in the previous sections. We will refer to this source 
of axion energy density as the misalignment energy density. The third source is the decay 
of cosmic axion strings [1S,2(|. In Ref. [19], it was found that the energy density resulting 



from this process is two orders of magnitude larger than the misalignment energy density, 
while in the estimate of Ref. ||2U|| , the energy density from the decay of cosmic strings 
is comparable to the misalignment energy density. At T ~ T\ the decay of strings 
will also produce an inhomogeneous axion field. While we can not describe the initial 
configuration emerging from string decay by the distribution of Eq. (|3.6|) , we expect that 
the attractive non-linear self interaction will also play a role here, and the evolution will 
proceed along the lines described in Sec. Ill and will result in high density peaks. The 
fourth potential source of axions is related to the collisions and subsequent disappearance 
of axion domain walls. In this section we discuss this last process. 

In most network of vacuum domain walls is a cosmological disaster, since 

they soon come to dominate the energy density of the Universe ||21|| . Fortunately, in the 



N = 1 axion model domain walls are effectively unstable and this problem is avoided. The 
process by which the domain walls disappear is through collisions of the string network 
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with the walls. The usual assumption is that when a string loop (with a wall on the 
inside) hits a large wall the two pieces of wall annihilate and a vacuum hole is produced 
in the large wall. Since the surface energy of the hole is smaller than the surface energy 
of the wall, the vacuum hole expands and devours the wall. For an infinite domain wall 
there will be roughly one hole in the wall per Hubble radius, so in a couple of Hubble 
times the holes quickly overlap and the wall disappears. Domain walls of finite size 
(size smaller than the Hubble radius) form closed surfaces and shrink by themselves. An 
oversimplified point of view would be that all of the energy released in the disappearance 
of the domain walls is transferred to relativistic axions, which subsequently redshift away 
and would become an insignificant source of axion energy density. 

We shall argue here that the hole in a domain wall formed by a string-loop intersection 
is not vacuum, but rather consists of a bound state of two pieces of domain wall (which 
in the simplified scenario annihilated each other) corresponding to a generalization of 
the "breather" solution of the (1 + l)-dimensional sine-Gordon model ||22|| . That is, the 
vacuum wall network is transformed into a breather wall network. The breather wall 
effectively evolves as a dust wall rather than a domain wall, so it represents cold dark 
matter. 

We consider here the axion field in Minkowski space-time with planar symmetry as 
a function of two coordinates, time Xq and one spatial direction, x\. It is convenient to 
introduce the dimensionless variables t = x§m a and z = X\m a . The relevant breather 
solution to the equation of motion 9 — Q" + sin 6 = (dot denotes d/dt and prime denotes 
d/dz) has the form | 22f 

sin(£ / r) 



9 v (t, z) = 4 arctan 



v cosh(,2 / L) 



(4.1) 



where r = y/1 + v 2 /v and L = yl + v 2 . One can interpret this solution as a bound state 
of two static domain walls (or kinks), #kink = ±4 arctan[exp(,z)]. The free parameter v 
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of the breather is related to the binding energy. Larger v corresponds to larger binding 
energy of the kinks. In terms of the spatial energy density distribution, the breather 
looks like a domain wall with an effective width L, but unlike the usual static domain 
wall, the field coherently oscillates with period r in the breather. When v — > the period 
tends to infinity, and when v — > oo the field oscillates with a frequency equal to the axion 
mass. The width of the breather scales in the opposite way with v: as v — > the width 
is m~ 1 , and the width grows proportional to v at large v. 

When considering wall like structures, it is convenient to introduce a surface stress- 
energy tensor of the wall, = /f^ T^ v dz. While all the components of in the 
breather solution oscillate with time, the surface energy density is constant: 



S° = 16/> a /VT+l/», (4.2) 

and S^ u is a diagonal tensor. The spatial components of S^ v are oscillating functions. 
However, when considering the macroscopic properties of a wall, the relevant quantities 
are averages over an oscillation period. Upon averaging over an oscillation period (S z z ) = 
0, as it must be for a wall of any nature PB[ . For the time-averaged surface tension we 
find 



(S\> = 8f 2 a m a 



2(v / TT 



v 2 — v 



V 



5% = SP. . (4.3) 



l + v 2 . 

As v — > 0, the stress-energy tensor tends to the vacuum stress tensor, with S° = S, 
where S° is twice the energy density of a single kink. However, at large v, we have 
S° ~ 16 /m a v and S ~ 8A^/m a v 3 . With increasing v, this tends to the stress- 
energy tensor of a dust shell. So in the expansion of the Universe, the surface density 
of the breather wall has to decrease and v has to grow. Using Eqs. ( |4.2| ) and ( |4.3j ), we 
obtain as a solution to the planar wall equations of motion p3| 



S°o(R) = 3V l m + a f > (4-4) 

20 



where we assumed constant m a and have have normalized the scale factor in such a 
way that R — 1 at the moment when v = 0. Note that the number of axions per 
unit area is conserved at large R. Despite the fact that the breather is a bound state, 
its surface energy density decreases in expansion, exactly as the energy density of the 
solution presented in Fig. 11. 

We can visualize the formation of the breather network in the following scenario. 
When domain walls form at T ~ T\, every string loop develops a wall inside (an "an- 
tiloop" develops a wall outside). When a string hits a large segment of wall, the in- 
tersection region will not be empty, but will be a bound state of two domain walls. In 
the idealized approximation of planar symmetry, the bound state will correspond to the 
breather solution of Eq. (}4.1|). However, since the perfectly planar situation is unrealistic, 



the question arises whether breather walls are stable. 

To answer this question we numerically integrated the axion equations of motion in 
Minkowsi space-time with initial conditions corresponding to a perturbed breather wall. 
We evolved an axisymmetric configuration 9 = 8(t,z,r), which initially corresponded 
to the breather field distribution of Eq. Q4.1| ) with v = v(r). The value of v, and the 
corresponding pressure, was larger in the center (note that Eq. (|4.3| ) corresponds to a 
system with negative pressure). In a sense, this configuration coresponds to a bubble of 
new phase of lower energy density, and it is expected to expand. The question is whether 
the field inside the "bubble" will tend to a breather solution with a new constant value of 
v as the boundary of the bubble propagates outward. We have found that this does not 
occur: the breather wall is unstable. However, the energy density in the breather does 
not dissipate, but the breather fragments into clumps very similar to those discussed in 
Sec. III.C. This result is not unexpected in view of the attractive nature of the axion 
self-interaction. The energy density profile in the r-direction is presented in Fig. 13 at 
several moments of time. 
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Our conclusion in this section is that the decay of the axion domain wall network can 
provide yet another channel for axiton production in the axion distribution. 

V. DISCUSSION 

In principle, all axion miniclusters could be relevant to laboratory axion search ex- 
periments, since even for $ as small as 1, the density is 10 10 times larger than the local 
galactic halo density [see Eq. ( |2.3|) 1. Moreover, as we have noted already, the energy 
density in an axiton after it separates out from the general expansion will be $ 4 times 
larger than the energy density at Teq. For example, a rather moderate density contrast 
of $ = 30 at Aqcd > T > Teq will correspond to roughly an additional factor of 10 6 in 
the energy density of the axiton at T <C Teq. 

The probability of a direct encounter with a minicluster is small. Let's assume that 
all of the axions end up in miniclusters of mass 10 _9 M Q , density 10 _14 g cm -3 , and radius 
4 x 10 12 cm. Using a local halo mass density of 5 x 10~ 25 g cm" 3 would give a minicluster 
number density of 7,000,000 pc~ 3 . With a typical velocity of 250 km s _1 the encounter 
rate would be 1 per 25 million years. Although the signal in an axion detector from a 
close encounter with a minicluster would be enormous, it might be a long wait. So the 
interesting question arises, could there be any other astrophysical consequences of very 
dense axion clumps? Below we shall discuss the possibility of "Bose star" formation in 
axion miniclusters. 

The physical radius of an axiton at Teq is larger by many orders of magnitude than the 
de Broglie wavelength of an axion in the corresponding gravitational well. Consequently, 
the gravitational collapse of the axion clump and subsequent virialization can be described 
in the usual terms of cold dark matter particles. In a few crossing times some equilibrium 
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distribution (presumably close to an isothermal distribution) of axions in phase space 
will be established. It is remarkable that in spite of the apparent smallness of axion 
quartic self-couplings, |A | = (f n / f a ) A ~ lCT 53 /^ 4 , the subsequent relaxation in an axion 
minicluster due to 2a — > 2a scattering can be significant as a consequence of the huge 
mean phase-space density of axions H- In the case of Bose-Einstein statistics the 



inverse relaxation time is (1 + n) times the classical expression, or r^ 1 ~ nv e ap a /m a , 
where a is the corresponding cross section. For particles bound in a gravitational well, 
it is convenient to rewrite this expression in the form |T2| 



rn 7 a \ a 2 p a 2 v 2 e . (5.1) 

The shallower the gravitational well for a given density of axions, the larger the mean 
phase space density, and consequently the smaller the relaxation time due to the v 2 
dependence in Eq. ( |5.1| ). Note also the dependence of the inverse relaxation time upon 
the square of the particle density. 

The relaxation time ( |5.1| ) is smaller then the present age of the Universe if the energy 
density in the minicluster satisfies 

Pio > 10V 8 /fo, (5.2) 

where pio = p/(10eV) 4 and Vs = f e /10~ 8 . If this occurs, then an even denser core 
in the center of the axion cloud should start to form. An analogous process is the so- 
called gravithermal instability caused by gravitational scattering. This was studied in 
detail for star clusters, where the "particles" obey classical Maxwell-Boltzmann statistics. 
Axions will obey Bose-Einstein statistics, with equilibrium phase-space density n(p) = 
n con d+[e /3 ' E — containing a sum of two contributions, a Bose condensate and a thermal 
distribution. The maximum energy density that non-condensed axions can saturate is 
Pthor ~ m t v ei "which corresponds to n t h er ~ 1- Consequently, given the initial condition 
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ri > 1, one expects that eventually the number of particles in the condensate will be 
comparable to the total number of particles in the region if relaxation is efficient. Under 
the influence of self gravity, a Bose star [23], [25[] then forms [|T2j] . One can consider a 



Bose star as coherent axion field in a gravitational well, generally with non-zero angular 
momentum |L|| . 

Comparing Eqs. ( |2.3|) and (|5.2|) , we conclude that the relaxation time is smaller than 
the present age of the Universe and conditions for Bose star formation can be reached in 
miniclusters with a density contrast $ > 30 at the QCD epoch. 

Under appropriate conditions stimulated decays of axions to two photons in a dense 



axion Bose star are possible [16, Q (see also [27]), which can lead to the formation of 
unique radio sources — axionic masers. In view of the results of this paper we conclude 
that the questions of axion Bose star formation, structure, and possible astrophysical 
signatures deserve detailed study. 

In conclusion, we have presented a 3-dimensional numerical study of the evolution of 
inhomogeneities in the axion field around the QCD epoch, including for the first time 
important non-linear effects. We found that the non-linear effects of the attractive self- 
interaction can lead to a much larger density of axions in miniclusters than previously 
estimated. Large amplitude density contrasts form solitons we call axitons, and resemble 
the bound-state "breather" solutions of the (1 + l)-dimensional sine-Gordon model. The 
increase in the axion density may be sufficiently large that axion miniclusters formed by 
the fluctuations might exceed the critical density necessary for them to relax to form 
Bose stars. 
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